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Abstract 

We consider a three dimensional system consisting of a large number of small 
spherical particles, which move due to gravity or with laminar shear and 
which merge when they cross. A size ratio criterion may be applied to re- 
strict merging to similar sized particles (locality of interactions) or particles 
dissimilar in size (nonlocality). We perform direct numerical simulations 
(DNS) of this particle system and study the resulting mass spectra. In mean 
field approximation, these systems can be described by the Smoluchowski co- 
agulation equation (SCE). DNS of the particle system with locahty enforced 
show the scaling solutions or Kolmogorov-Zakharov spectra for the SCE, sig- 
nifying a constant mass flux. DNS without a size ratio criterion show —4/3 
scaling for large particles in a system with gravity, signifying a constant flux 
in number of particles, which we also find analytically by assuming nonlo- 
cality of interactions in the SCE. For laminar shear, this nonlocality is only 
marginal, and our DNS show that a correction to the scaling solution is 
required. 

1 Introduction 

Systems of coagulation are abundant in the physical sciences, with notable 
examples in the atmosphere [HE], and applications in chemistry [3111] and 
astrophysics [5] . Our initial interest in these systems was raised by their use 
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in the theory of warm rain formation [HI El IE] and the existence of scahng 
solutions [31 E] • 

We consider coagulation as an event where two particles of mass mi and 
m2 merge to form a particle of mass m = mi + m2. The exact nature 
of how these particles meet and merge is described by a collision kernel 
K{mi,m2)- This description is reminiscent of the turbulence formalism, 
where wave interaction between wave numbers ki and k2 is described by an 
interaction coefficient Vki2 [ID]. We wish to continue this analogy, and create 
a particle system with a turbulence setup, i.e. with particles forced at some 
small mass mo, and removed once they reach a large mass sink M. Assuming 
this system reaches a steady state, we hope to find scaling solutions conform 
Kolmogorov-Zakharov solutions in turbulence. 

From a mathematical perspective, particle coagulation has been discussed 
with a focus on scaling and self-preserving solutions to the Smoluchowski 
coagulation equation (SCE), which was originally derived as a mean field 
approximation to the physical coagulation processes for Brownian motion and 
laminar shear [TT] . Scaling solutions for general kernels of coagulation have 
been classified, with a particular interest for unphysical solutions that imply 
gelation, i.e. the formation of a particle of infinite mass at finite time [T2l[T3l 
[H]. Similarities between the SCE and turbulence theory have been discussed 
previously [15] , with an explicit discussion of locality of interactions. In this 
paper, we use results from an interacting particle system to highlight the 
importance of this assumption of locality when deriving scaling solutions. 

First, in Sect. [2] we will introduce the interacting particle system, which 
does not assume mean field behaviour or any other properties of the SCE. 
The system can model particles moving due to gravity, a coagulation process 
known as differential sedimentation [3], or particles moving due to laminar 
shear. Using a simple size ratio criterion, we can ensure only similar-sized 
particles merge, making sure the locality assumption is met, or we ignore the 
size ratio and allow all particle interactions. We also perform simulations in 
which we only allow interactions between particles very dissimilar in size to 
study nonlocal interactions more directly. 

We continue by introducing the SCE and provide the well-known scaling 
results (Sect. [3]). We comment on further properties of the kernels for laminar 
shear and differential sedimentation, in particular the convergence of the 
collision integral in the SCE (Sect. H]). 

In Sect. [5] we revisit the approximate solution for the system with dom- 
inance of nonlocal interactions in differential sedimentation [T^]. Using the 
IPS, we retrieve the scaling solution. We find that the system without forced 
locality has a bottleneck and that it shows non-mean field behaviour. We 
provide several techniques to reduce this finite size effect and compare results 
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to the nonlocal solution. 

Finally, we study laminar shear in Sect. El which in terms of the SCE is 
known to be a boundary between gelling and non-gelling systems. We do not 
provide a nonlocal solution for this case, but we present numerics from the 
IPS showing a deviation from the scaling solution that clearly is a result of 
nonlocal interactions. 

2 The interacting particle system 

We have developed a very simple system of a large number of interacting 
particles (IPS) with which we can easily test a variety of properties of coag- 
ulation. We use a setup typical for turbulence simulations, in which we force 
mass at small scales and dissipate at large scales, and let the mass or volume 
spectrum develop. Q Mass transfer occurs because of particle coagulation, 
where two particles of volumes mi and m2 form a new particle of volume 
m = nil + "^2- We may allow or disallow a coagulation event between two 
particles depending on their size ratio mi/m2, thus enforcing either locality 
or nonlocality of interactions. If we allow all particle interactions, we refer 
to this system as "free merging" . 

In the IPS, each time step we introduce a fixed number A^o of particles 
of discrete multiples j of the minimum particle volume tuq with j e [1,5], 
uniformly distributed in space. This space is a periodic box of size x 
Ly X Lz, where typically = Ly = = 10cm. We consider the particles 
spherical with minimum radius tq = 0.01cm and to have no interaction with 
the surrounding "fluid". Particle coagulation in the IPS happens when the 
distance between the centers of two particles of volume mi and m2 is smaller 
than the sum of their radii r = ri + r2. A new particle of volume m = mi+m2 
is then created at the center of mass of the original particles. If the rules of 
collision disallow this event, the particles will cross one another, having no 
effect on each other's trajectory. Finally, we typically choose a strict cutoff 
M = lO'^mo, so that any particles formed that have m > M are removed 
from the system at the end of the time step. 

In a system with gravity, large particles will overtake smaller ones because 
of their higher terminal velocity, and thus coagulation may occur. This 
system is referred to as differential sedimentation We assume particles 
move at Stokes terminal velocity u ~ Gr^ ~ m^/^, where we choose G = 2 ■ 
10^cm~^s~^, so that a particle of radius 0.1cm travels at 20cm -s"^, equivalent 
to an air bubble rising in water [TT] . 

"'^We will make no formal distinction between mass and volume in this paper. 
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For laminar shear, we choose a shear function f{z) = jz with 7 = 12.5s~^, 
so that those particles at z = 10 travel at 125cms~^. Particles ai z' = z + dz 
will overtake those at z, and coagulation may occur. To ensure periodicity 
between particles near z = and those near z = Lz, when we check for these 
events we move all particles with < z' < 0.5 to z' + Lz- The x-position of 
those particles is then altered so that Xnew = ^om + iTLz, where T is the time 
that has elapsed since the start of the simulation. Once we have checked for 
collisions over this periodic boundary, we move the particles back to z' and 
transform Xnew back to Xoid- 

Simulations for a system with free merging will have no collision rule ap- 
plied. In systems with forced locality or nonlocality, we check the size ratio 
mi/ 1712 between two overlapping particles. If this satisfies 1/q < mi/m2 < q, 
where typically g = 2, we have locality. If we only allow particle coagulation 
when the above criterion is satisfied, we call this a system with "forced lo- 
cality" . If we only allow particle coagulation when the above criterion is not 
satisfied, we call this a system with "forced nonlocality" . 

At long times, we assume the IPS has reached a statistical steady state, 
and we produce the particle size spectrum. We divide the volume range 
into logarithmic bins bk with G N such that bk = [l.l'^mo, l.l'^'^^mo). At 
each time step, we add all particles to their corresponding bin bk and after 
10, 000 time steps we average to find the particle size spectrum. This method 
smooths the resulting spectrum and has no noticeable effect on the scaling. 

3 Mean field approximation 

The Smoluchowski coagulation equation [TT] has been applied to numerous 
aggregation and clustering mechanisms in physics and chemistry [H [3l [181 E] • 
We will check whether it is an appropriate mean field approximation for the 
IPS. We write the SCE 



where n{m) is the density of particles of volume m. The function K{mi, m2) 
is the collision kernel for particles mi and m2, which describes the collision 
process between particles merging to form a particle of size mi + m2, and u 
is the particle velocity. 

In past research, scaling solutions to the SCE have been derived using 
dimensional analysis [31 These solutions have also been derived by ap- 
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plication of the Zakharov transformation to the coUision integral \i9l, and 
conform turbulence terminology they can be interpreted as a Kolmogorov- 
Zakharov cascade of volume or mass [20l E]. For differential sedimentation 
this solution is 

n{m) ~ m"^^/^ , (2) 

and for laminar shear 

n{m) ~ . (3) 

Both the scaling argument and the Zakharov transformation rely on locality 
of interactions, i.e. convergence of the collision integral in the SCE. However, 
we will now show that both the kernel for differential sedimentation and the 
one for laminar shear have scaling exponents that lead to divergence. 



4 Convergence of the collision integral 

Let us rewrite the SCE ([T]) for a system with a small and large mass cutoff, 
i.e. a truncated spectrum |21j 



— — h u ■ Vn = / dmiK(mi,m — mi)n(mi)n(rn — rrii) (4) 

dt Jmo 

pM 

— / dmiK(m,mi)n(m)n{mi) , 

Jmo 

where we have used the symmetry of the first integral. 

We assume a scaling solution to this equation of the form n{m) ~ m~^, 
which requires locality of interactions, or collision of the integral as rriQ | 
and M oo. Also, we assume two types of scaling behaviour of the kernel, 
namely 

K{mxi,mx2) = ■m^K{xi, X2) , (5) 

and for xi ^ X2 

K{xi,X2) r^x^x'^ . (6) 
Let us introduce m_ and m_|_ such that 

rriQ <C m_ <C m <C M , (7) 
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and we rewrite eq. using these new delimiters 




/C(mi, m — nil) — /C(mi, m)dmi 
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I'm/2 

+ / /C(mi, m — mi) — /C(mi, m)dmi 



J m_ 



f /C(mi,m)dmi 

m/2 

f /C(mi, m)dmi , 



where we introduce /C(mi,m2) = K{mi^m2)n{mi)n{'m2) and combine the 
two integrals where possible. 

Clearly, if we introduce an additional kernel £{mi, 1712) to force locality in 
the system we can ignore the two integrals in eq. ([H]) that involve mo and M, 
and the collision integral is convergent. That is, when £^ = if mi/m2 < l/o' 
or mi/m2 > q with g > 1 finite, then eventually mo/m < 1/q and M/m > q 
so taking the limits mo i and M ^ 00 will not lead to divergence. 

We leave m_, m, m+ fixed, and study the limit mo J, 0. Because of the 
inequality ([7]), we can use a Taylor expansion on the first integral in eq. ([8]) 
and write 



Thus, for the integral (l9j) to converge as mo J, 0, we require x < 2. 

Under the same conditions, we can write the last integral in eq. (IHl) as 



and if we let M — ^ 00 we require x > A + 1 for convergence of the integral. 
Both requirements x < 2 and x > A + 1 have been derived previously [20] ■ 

Now, for differential sedimentation we have A = 4/3 and x = 13/6, which 
implies divergence of the collision integral, and we require the additional 
kernel S with finite q to force locality in this system so that we can achieve 
the scaling solution. For laminar shear we have A = 1 and x = 2, which puts 
it on the boundary of convergence. 

The kernel of differential sedimentation is traditionally regarded as un- 
physical, because the growth rate of large particles increases faster than their 
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actual volume [12]. It has been shown that by a simple modification of 
this kernel, scaling solutions can be realised in numerics with a particle sys- 
tem [16]. Both kernels, or at least their exponents A have also been discussed 
in relation to gelation, i.e. the formation of an infinite mass particle at finite 
time [201 1121 1131 E] • Within this family of kernels, the exponent A = 1 forms 
a boundary between gelling and non-gelling systems [221 123]. The introduc- 
tion of a finite large size cutoff limits the possibility of gelation in the IPS. 
However, the effect of a finite cutoff on the process of gelation in the SCE is 
unclear. 



5 Differential sedimentation 



Numerical simulations of an interacting particle system have shown that a 
free merging system for differential sedimentation does not show the pre- 
dicted scaling eq. ([2]) [H]. Instead, the particle density spectrum retrieved 
from the simulations can be explained by rewriting the SCE and assuming 
nonlocality of interactions to be the dominant process (see Appendix |X] and 
Ref. [IE]). A similar approach was taken for the study of nonlocality in drift 
waves in turbulence [21] and the evolution of galaxies [25] . 

For differential sedimentation, we may write this approximate nonlocal 
solution to the SCE ([1]) as 



n m 



no 



m 



mo 



-4/3 



exp 



WO) I ( ^ 



-1/3 



'11^ 



with uq = n(mo), and W{Q) = W (^2^ j the Lambert W function or product 

log. Since we assume M ^ mo, we may approximate W{Q) by logB. 

The dependence of the solution ( ITIl) on sharp cutoffs at mo and M and 
the traditional solution of "instantaneous gelation" [T2] of similar kernels 
requires a word of caution concerning its validity as a steady state solution 
to the SCE. We will however check its applicability as an approximation to 
the IPS with cutoff. 

Note that at large m, the solution fllip follows the scaling n{m) ~ m~^/^. 
This limit in itself is a conservation of number of particles for large m, i.e. 
there is a cascade for the number of large particles in the system. Since we 
do not use a strict monomer source mo we will not focus on a fit of eq. (fTTj) 
towards smaller m. 

We have performed numerical simulations for a system with differential 
sedimentation under different rules of collisions, as described in Sect. [21 The 
data for forced locality in Fig. [1] shows rather strong fluctuations within the 



7 




Figure 1: Particle size density spectra for differential sedimentation. Circles: 
forced locality; upward triangles: free merging; downward triangles: forced 
nonlocality. Dotted lines show the —13/6 and —4/3 power laws. Dashed 
line signifies the end of the forcing range at 5mo; cutoff for all simulations is 
at M = lO^mo. In all figures of mass spectrum, markers only indicate each 
graph, not all individual data points. 
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Figure 2: Positions of all particles at a given time step in a vertical (x, z)- 
slice of width 0.5cm of the domain for a free merging system with differential 
sedimentation and cutoff M = lO^mo- Small dots are particles of size 1 — 
25mo; intermediate are of size 25 — 625mo; large are of size larger than 625mo. 



9 



Figure 3: Similar setup as Fig. [2] but with random reintroduction at 2; = 0. 



cascade, which we blame on the actual mechanism of differential sedimenta- 
tion, in which particles of the same size will not merge. Thus, if we force 
particles of integer multiples j of mg uniformly in the range j G [1, 5], we have 
particles formed of size 3 and 5 due to coagulation in addition to the forcing, 
and we expect peaks in their distribution. Similarly, particles of size 4 coag- 
ulate with sizes 2, 3, and 5, but are not created (apart from forcing), and we 
expect a below average value for their number. These initial inhomogeneities 
are consequently carried through the cascade. 

Also in Fig. [T] we plot the data for forced nonlocality and free merging. 
We clearly see that the free merging and nonlocal spectra do not follow a 
strict scaling law. Instead, we note that the simulation for free merging 
shows a spectrum similar to that of forced nonlocality, but neither can be 
convincingly fitted to a —4/3 power law at the large m as predicted by the 
solution ( fTTj) . We see an overestimation of number of large particles closer 
to a spectrum n(m) ~ m~^, which indicates a finite size effect or bottleneck 
in our simulation. 

A snapshot of the free merging system at large times shows that large 
particles have dug out channels void of any particles (Fig. [2]). Due to peri- 
odicity of the system, the large particle will have traversed the whole length 
in a time too short for smaller sized particles to hang around and not be 



10 



merged with the larger particle. This causes stagnation of the mass flux and 
consequently a bottleneck accumulation on particle mass distribution. As 
we see, such a bottleneck mechanism has a non-mean field origin. We try 
to avoid this behaviour, which contradicts the mean field assumption used 
in the derivation of the SCE, by stretching the domain in the direction of 
motion. The result for such a system is shown in Fig. HI We note that 
the bottleneck effect is not as strong as in the cubic domain, and the large 
particle behaviour is close to the nonlocal solution ( fTTll . 

Independently, we may try to avoid the non-mean field behaviour by 
adapting the periodicity of the system, effectively introducing mixing at a 
particular level of z. When a particle leaves the system at -|- z\ where 
Lz is the length of the domain in the direction of motion, we reintroduce 
the particle at z\ but with new and random x, ^/-coordinates. The resulting 
spectrum for this system is also shown in Fig. HI and a snapshot is given 
in Fig. [31 We note that the channels for large particles are not throughout 
the system anymore, and that there are typically fewer large particles. In- 
terestingly, the spectrum is not significantly different from the spectrum for 
the stretched domain. A similar modification of periodicity has been used 
previously for an IPS with laminar shear [26] . 

In an attempt to gain a longer range of n(m) ~ m~^/^ for large m, we 
compare the data for the IPS with data from direct numerical simulations 
of the SCE with cutoff ([1|). The results from Fig. [3] have been realised 
with forcing of equal volume into both ttiq and Im^. A wider inertial range 
has been realised with the SCE numerics, and we see the same qualitative 
behaviour of the spectrum with —4/3 at large m for all simulations. 

From Fig. [5] we see good agreement between the analytic solution ( [TTl) 
and the SCE numerics. This is not entirely surprising, since this solution was 
derived from the SCE with cutoff, but using the assumption of dominance 
of nonlocal interactions. Both the SCE numerics and the analytic solution 
overestimate the IPS of which they are an approximation. The deviation 
of the SCE from the IPS may be due to non-mean field behaviour and the 
finite domain, or it could be related to gelation. The exact nature of this 
deviation will require further study of the IPS and the influence of the source 
and cutoff. 

6 Laminar shear 

For laminar shear we cannot find a nonlocal solution of the form f[TTl) (see 
Appendix [A]). Instead, we present the results from the IPS and look for 
differences between results for the different rules of collision. 
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Figure 4: Free merging simulations with cutoff M = lO^^mo- Circles: cubic 
domain; upward triangles: cubic with random reintroduction; downward tri- 
angles: no random reintroduction, but domain size 5 x 5 x 40. Dotted lines 
show the —13/6 and —4/3 power laws. Dashed hne signifies the end of the 
forcing range at 5mo. 
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Figure 5: Free merging simulations with increasing cutoff M. M = lO^mo 
in cubic domain and with random reintroduction (circles). Solid line: SCE 
numerics with M = lO^mo- Dashed line: nonlocal solution ffTTl) . Inset: SCE 
numerics and nonlocal solution for M = lO^^mo- Spectra shown are for 
n{m) = n{m)/n{mo)- Dotted line shows the —4/3 power law. SCE numerics 
courtesy Colm Connaughton. 
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The particle spectra for all three different setups are shown in Fig. El 
and we note that the simulation for forced locality clearly follows the scaling 
solution n{m) ~ m~^. Also, the fluctuations down the spectrum that were 
so clear for differential sedimentation are not apparent here, so we assume 
that these were indeed due to the specific kernel. 




Figure 6: Particle size density spectra for laminar shear. Circles: forced 
locality; upward triangles: free merging; downward triangles: forced nonlo- 
cality. Dotted line shows the —2 scaling solution. Dashed line signifies the 
end of the forcing range at 5mo; cutoff for all simulations is at M = lO'^mo. 



Again, the free merging system is similar to the nonlocal system in the 
particle spectrum, and we note that near the source and sink they deviate 
from the —2 scaling. The absence of this deviation in the data for forced 
locality clearly indicates that this is a result of nonlocality. In Fig. [7| we 
show the compensated spectra m^n(m) for all simulations. We see that the 
deviation from the —2 power law for free merging and forced nonlocality is 
a slowly evolving function of mass. 

The results from Fig. [7| combined with the knowledge that the exponent 
A = 1 for the SCE is critical in terms of convergence and gelation leads 
us to believe we should look for corrections to the —2 spectrum. Following 
Kraichnan's work in two-dimensional turbulence [271 EH], where he predicted 
logarithmic corrections for the enstrophy cascade which also has marginal lo- 
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Figure 7: Compensated particle size density spectra, rn?n{m) for the same 
setup as Fig. [6l 



cality, we assume the corrections to the —2 scaling solutions to be logarithmic 
as well. 



7 Discussion 

We have performed numerical simulations of an interacting particle system 
(IPS) under different sets of rules for both differential sedimentation and lam- 
inar shear. The system does not assume mean field behaviour, yet provides 
results that are remarkably consistent with theoretical predictions using the 
SCE for the respective mechanisms of coagulation. 

The IPS with forced locality confirms the KZ scaling laws for the SCE, 
whilst a free merging system deviates from said solutions, and its spectrum 
is more similar to a system that assumes nonlocal interactions to be domi- 
nant. For differential sedimentation we have predicted a solution which has 
a particle number cascade of —4/3 at large m. Our numerics show spectra 
close to this behaviour, but the simulations suffer from a finite size effect. 
In the standard system, due to periodicity large particles dig out channels 
which become devoid of smaller particles, which slows down the mass cas- 
cade. This in turn leads to a bottleneck accumulation of mass at large m. 
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which behaviour is non-mean field in origin. We hmit this effect by stretching 
the system or by introducing mixing, both of which numerics show a —4/3 
spectrum for large m. 

The exponent A = 1 in the kernel for laminar shear is known to be 
marginal in terms of convergence of the collision integral f20] and we recog- 
nise this marginal nonlocal behaviour in the results from our IPS. In a com- 
pensated spectrum we note that corrections to the —2 scahng are necessary, 
and we will investigate such corrections in a future paper. 

Traditionally, both differential sedimentation and laminar shear are mech- 
anism thought of as leading to "instantaneous gelation" fi2[ [13] . We note 
that in spite of this prediction of unphysical behaviour, we find a persistent 
solution with our IPS, which can be described by an approximate solution 
using the SCE with cutoffs. 

Finally, we note that the remarkably good agreement of the IPS with the 
predictions for the SCE with cutoffs has also been confirmed by simulations 
done for self-similarity ^J, and the nonlocal prediction of n(m) ~ m~'^^^ for 
large mass can also be recognised in SCE numerics. 
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A Nonlocal solutions to the collision integral 

We continue the analysis in Sec. H] from eq. ([HI), forgetting the assumption 
of a scaling solution for n(m), and instead assume that nonlocal interactions 
are dominant in the system. Again leaving m_, m, m+ fixed, we can use 
eq. ([9]) and eq. (ITOl) to write 



where we replaced by n{m), and we have replaced the integrals in eq. (jH]) 
and ( ITOl) by ci and C2 respectively, again replacing by n{m). 

Now, we still assume that the inequahties (JTj) hold, but we keep uiq and 
M fixed. Assuming this system may reach steady state, we obtain an ODE 
for eq. ffT^ which for A > 1 has the solution 



Dn 



Ci(9m(m n{m)) — C2n{m) 



(12) 



Dt 



n{m) = Cm exp 



1 C2 

A - Ici 



l-A 



(13) 
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We will not use this approach to find a solution for the marginal case with 
A = 1. 

If we define no = n{mo), and set m = mo in eq. ( fT3l) . we may write 



C = TioTriQ exp 
and we write for the solution ( JT3l) @ 



1 C2 1 



1 - Aci 



rrin 



(14) 



nim) = no — exp 
V m / 



1 £i, 

A - Ici 



m 



(15) 



We may also find approximations for C2 and Ci in order to find a steady state 
solution for A > 1. Using eq. (fT5|) we write 



/M 



A _1-A 



niQ ")] dm 



1 ) 



(16) 



with A = 1^^- Since we have M, m ^ mo, and 1 — A < 0, we may write 

C2 nom^Mexp [-v4mo~^] . (17) 

For ci we write 

ci = / nom^ml^^exp [A(m}"^ - mg"^)] dmi , (18) 

J mo 

and we note that as mi — * mo, the term exp[...] —>■ 1. Thus, we may write 



ci f» nomo 



2- A 



-m 



,2-A 



2- A 



notriQ 



(19) 



where we note the additional restriction that A < 2. 
Now we can write 



1 C2 

A - Ici 



2 - A 



A = mo~^Mexp [-Aml"^] 

A — 1 



(20) 



Let us write A = Om^ ^, and find 



, 2-XM 

e = — e 

A — 1 mo 



(21) 



^Alternatively, one can match the constant C to the flux of particles near mo- 
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or 

which has as a solution for 9 

e = w 



2-XM 



A — 1 mo ' 
2-XM 



A — 1 mo 



W{Q) , 



(22) 



(23) 



where W{z) is the Lambert W function or product log. This is a single- 
valued, real function under our assumption 1 < A < 2. 

We find our general solution to the Smoluchowski equation with nonlo- 
cality and homogeneity A > 1 to be 



n(m) ~ nomQiri ^ exp 



w{e) I ( ^ 

mo 



l-A 



(24) 



Note that such a solution should be taken with caution, as it uses the assump- 
tion of sharp cutoffs mo and M. As such, the exponential "spike" towards 
small m will depend on the forcing near mo- We do expect the large m 
scaling of n{m) ~ m"^ to be a verifiable feature of the nonlocal solution. 
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